################################################
## CONJOINT ANALYSIS - NEW ZEALAND  ###
################################################


rm(list=ls(all=TRUE))
detach(US.conj)


### Load libraries
library(arm)
library(stargazer)
library(foreign)
library(sjPlot)
library(ggplot2)
library(ggthemes)
library(ggeffects)
library(foreign)
library(mediation)
library(cregg)
library(FindIt)


### Read and attach data

#Find your path

NZ.conj <- read.csv("C:/Users/gmagni/OneDrive/SideProjects/3-AndyNEW/00-Outlin&Paper-Appendix/00-Latest/00-JOP/R&R/FinalSubmission/Data/NZ-LGBT-Conj.csv")


attach(NZ.conj)
n = nrow(NZ.conj)




#################################################




### Code variables


## Conjoint variables: candidate choice

table(CJ.1.Rat.1_1)
rating.B.nz <- as.numeric(as.character(CJ.1.Rat.1_1))
table(rating.B.nz)

table(CJ.1.ChoiceDummy)
dummy.choice.B.nz <- as.numeric(as.character(CJ.1.ChoiceDummy))
table(dummy.choice.B.nz)

table(CJ.1.REC_1)
race.B.nz = rep(NA,n)
race.B.nz[CJ.1.REC_1=="White"]="White"
race.B.nz[CJ.1.REC_1=="Maori-Pacific Islander"]="Maori-Pacific"
table(race.B.nz)

race.B.nz.factor <- factor(race.B.nz, levels = c("White", "Maori-Pacific"))
table(race.B.nz.factor)


table(CJ.1.REC_2)
sexorient.B.nz = rep(NA,n)
sexorient.B.nz[CJ.1.REC_2=="Straight"]="Straight"
sexorient.B.nz[CJ.1.REC_2=="Gay"]="Gay"
table(sexorient.B.nz)

sexorient.B.nz.factor <- factor(sexorient.B.nz, levels = c("Straight", "Gay"))
table(sexorient.B.nz.factor)


table(CJ.1.REC_3)
religion.B.nz = rep(NA,n)
religion.B.nz[CJ.1.REC_3=="Not religious"]="Not religious"
religion.B.nz[CJ.1.REC_3=="Christian"]="Christian"
religion.B.nz[CJ.1.REC_3=="Muslim"]="Muslim"
table(religion.B.nz)

religion.B.nz.factor <- factor(religion.B.nz, levels = c("Not religious", "Christian", 
                                                         "Muslim"))
table(religion.B.nz.factor)


table(CJ.1.REC_4)
age.B.nz = rep(NA,n)
age.B.nz[CJ.1.REC_4=="35"]="35"
age.B.nz[CJ.1.REC_4=="44"]="44"
age.B.nz[CJ.1.REC_4=="56"]="56"
age.B.nz[CJ.1.REC_4=="71"]="71"
table(age.B.nz)

age.B.nz.factor <- factor(age.B.nz, levels = c("35", "44", "56", "71"))
table(age.B.nz.factor)


table(CJ.1.REC_5)
health.B.nz = rep(NA,n)
health.B.nz[CJ.1.REC_5=="In a wheelchair since birth"]="Weelchair since birth"
health.B.nz[CJ.1.REC_5=="Healthy"]="Healthy"
health.B.nz[CJ.1.REC_5=="Overweight, has diabetes"]="Overweight, diabetes"
health.B.nz[CJ.1.REC_5=="HIV positive"]="HIV+"
health.B.nz[CJ.1.REC_5=="HIV positive since birth"]="HIV+ since birth"
table(health.B.nz)

health.B.nz.factor <- factor(health.B.nz, levels = c("Healthy", "HIV+", "Weelchair since birth", 
                                                     "Overweight, diabetes",
                                                     "HIV+ since birth"))
table(health.B.nz.factor)


table(CJ.1.REC_6)
polexp.B.nz = rep(NA,n)
polexp.B.nz[CJ.1.REC_6=="No previous experience"]="No experience"
polexp.B.nz[CJ.1.REC_6=="Town council member"]="Town council"
polexp.B.nz[CJ.1.REC_6=="Member of the House of Representatives"]="House of Repres"
table(polexp.B.nz)

polexp.B.nz.factor <- factor(polexp.B.nz, levels = c("No experience", 
                                                     "Town council", 
                                                     "House of Repres"))
table(polexp.B.nz.factor)


table(CJ.1.REC_7)
educ.B.nz = rep(NA,n)
educ.B.nz[CJ.1.REC_7=="Less than high school"]="No high school"
educ.B.nz[CJ.1.REC_7=="High school degree"]="High school"
educ.B.nz[CJ.1.REC_7=="Bachelor's degree"]="Bachelor"
educ.B.nz[CJ.1.REC_7=="Master degree"]="Master"
table(educ.B.nz)

educ.B.nz.factor <- factor(educ.B.nz, levels = c("No high school", "High school", 
                                                 "Bachelor", "Master"))
table(educ.B.nz.factor)


table(CJ.1.REC_8)
gender.B.nz = rep(NA,n)
gender.B.nz[CJ.1.REC_8=="Man"]="Man"
gender.B.nz[CJ.1.REC_8=="Woman"]="Woman"
gender.B.nz[CJ.1.REC_8=="Transgender"]="Transgender"
table(gender.B.nz)

gender.B.nz.factor <- factor(gender.B.nz, levels = c("Man", "Transgender", "Woman"))
table(gender.B.nz.factor)


## Causal mechanisms

table(CJ.B1.CM_1.Dummy) #left-leaning
cm.left.leaning = rep(NA,n)
cm.left.leaning[CJ.B1.CM_1.Dummy=="0"]=0
cm.left.leaning[CJ.B1.CM_1.Dummy=="1"]=1
table(cm.left.leaning)

table(CJ.CM_5.Dummy) #prefer as neighbor
cm.pref.neighb = rep(NA,n)
cm.pref.neighb[CJ.CM_5.Dummy=="0"]=0
cm.pref.neighb[CJ.CM_5.Dummy=="1"]=1
table(cm.pref.neighb)

table(CJ.CM_6.Dummy) #better chances to win election
cm.electorab = rep(NA,n)
cm.electorab[CJ.CM_6.Dummy=="0"]=0
cm.electorab[CJ.CM_6.Dummy=="1"]=1
table(cm.electorab)


## Controls

table(IncomeH)
income.h.nz = rep(NA,n)
income.h.nz[IncomeH=="No income"]=0
income.h.nz[IncomeH=="Less than  $5.000"]=1
income.h.nz[IncomeH=="$5.000-$10.000"]=2
income.h.nz[IncomeH=="$10.000-$15.000"]=3
income.h.nz[IncomeH=="$15.000-$20.000"]=4
income.h.nz[IncomeH=="$20.000-$25.000"]=5
income.h.nz[IncomeH=="$25.000-$30.000"]=6
income.h.nz[IncomeH=="$30.000-$40.000"]=7
income.h.nz[IncomeH=="$40.000-$50.000"]=8
income.h.nz[IncomeH=="$50.000-$60.000"]=9
income.h.nz[IncomeH=="$60.000-$70.000"]=10
income.h.nz[IncomeH=="$70.000-$85.000"]=11
income.h.nz[IncomeH=="$85.000-$100.000"]=12
income.h.nz[IncomeH=="$100.000-$150.000"]=13
income.h.nz[IncomeH=="$150.000-$200.000"]=14
income.h.nz[IncomeH=="$200.000-$250.000"]=15
income.h.nz[IncomeH=="More than $250.000"]=16
table(income.h.nz)

table(Educ)
educ.nz = rep(NA,n)
educ.nz[Educ=="No high school"]=0
educ.nz[Educ=="Some high school"]=1
educ.nz[Educ=="High school graduate"]=2
educ.nz[Educ=="ITP certification"]=3
educ.nz[Educ=="Some university but no degree"]=4
educ.nz[Educ=="Bachelor's or Honours degree"]=5
educ.nz[Educ=="Postgraduate degree (for example Master, PhD, ...)"]=6
table(educ.nz)

table(GenderID)
female.nz = rep(NA,n)
female.nz[GenderID=="Woman"]=1
female.nz[GenderID=="Man"]=0
table(female.nz)

table(YOB)
age.nz<-2018-YOB
table(age.nz)


table(PolIdeol.Eco_1)
eco.right.nz <- as.numeric(as.character(PolIdeol.Eco_1))
table(eco.right.nz)
prop.table(table(eco.right.nz))

table(PolIdeol.Soc_1)
soc.conserv.nz <- as.numeric(as.character(PolIdeol.Soc_1))
table(soc.conserv.nz)
prop.table(table(soc.conserv.nz))

table(Vote.Y.N)
voted.yes = rep(NA,n)
voted.yes[Vote.Y.N=="No, I didn't vote"]=0
voted.yes[Vote.Y.N=="Yes, I voted"]=1
table(voted.yes)

table(VoteChoice)
party.choice.nz = rep(NA,n)
party.choice.nz[VoteChoice=="ACT"]="ACT"
party.choice.nz[VoteChoice=="National"]="National"
party.choice.nz[VoteChoice=="Green"]="Green"
party.choice.nz[VoteChoice=="Labour"]="Labour"
party.choice.nz[VoteChoice=="NZ First"]="NZ First"
party.choice.nz[VoteChoice=="Other"]="Other"
#party.choice.nz[VoteChoice==""]="No Vote"
table(party.choice.nz)
prop.table(table(party.choice.nz))


table(KnowLGBT)
know.lgbt = rep(NA,n)
know.lgbt[KnowLGBT=="Yes"]=1
know.lgbt[KnowLGBT=="No"]=0
table(know.lgbt)

table(SexOrient)
lgbt.self = rep(NA,n)
lgbt.self[SexOrient=="Gay" | SexOrient=="Lesbian" | SexOrient=="Bisexual" | SexOrient=="Other"]=1
lgbt.self[SexOrient=="Straight"]=0
table(lgbt.self)

table(ReligY.N)
religious.yes = rep(NA,n)
religious.yes[ReligY.N=="Yes"]=1
religious.yes[ReligY.N=="No"]=0
table(religious.yes)

table(ReligID)
relig.ID = rep(NA,n)
relig.ID[ReligID=="Anglican"]="Anglican"
relig.ID[ReligID=="Presbyterian, Congregational, or Reformed"]="Presbyterian"
relig.ID[ReligID=="Hindu"]="Hindu"
relig.ID[ReligID=="Other Christian"]="O. Christian"
relig.ID[ReligID=="Catholic"]="Catholic"
relig.ID[ReligID=="Muslim"]="Muslim"
relig.ID[ReligID=="Methodist"]="Methodist"
relig.ID[ReligID=="Buddhist"]="Buddhist"
relig.ID[ReligID=="Other"]="Other"
table(relig.ID)

table(Religiosity)
religiosity = rep(NA,n)
religiosity[Religiosity=="Never"]=0
religiosity[Religiosity=="Only on special holidays"]=1
religiosity[Religiosity=="At least once a month"]=2
religiosity[Religiosity=="Once a week"]=3
religiosity[Religiosity=="Every day"]=4
table(religiosity)




## Add coded variables to dataset


NZ.conj$soc.conserv.nz <- soc.conserv.nz
NZ.conj$job.insec.nz <- job.insec.nz
NZ.conj$income.h.nz <- income.h.nz
NZ.conj$educ.nz <- educ.nz
NZ.conj$female.nz <- female.nz
NZ.conj$age.nz <- age.nz
NZ.conj$party.choice.nz <- party.choice.nz
NZ.conj$unemployed.nz <- unemployed.nz
NZ.conj$att.check.passed <- att.check.passed
NZ.conj$party.choice.nz <- party.choice.nz
NZ.conj$relig.ID <- relig.ID
NZ.conj$age.nz <- age.nz

NZ.conj$rating.B.nz <- rating.B.nz
NZ.conj$dummy.choice.B.nz <- dummy.choice.B.nz
NZ.conj$age.B.nz <- age.B.nz
NZ.conj$age.B.nz.factor <- age.B.nz.factor
NZ.conj$religion.B.nz <- religion.B.nz
NZ.conj$religion.B.nz.factor <- religion.B.nz.factor
NZ.conj$polexp.B.nz <- polexp.B.nz
NZ.conj$polexp.B.nz.factor <- polexp.B.nz.factor
NZ.conj$health.B.nz <- health.B.nz
NZ.conj$health.B.nz.factor <- health.B.nz.factor
NZ.conj$sexorient.B.nz <- sexorient.B.nz
NZ.conj$sexorient.B.nz.factor <- sexorient.B.nz.factor
NZ.conj$educ.B.nz <- educ.B.nz
NZ.conj$educ.B.nz.factor <- educ.B.nz.factor
NZ.conj$race.B.nz <- race.B.nz
NZ.conj$race.B.nz.factor <- race.B.nz.factor
NZ.conj$gender.B.nz <- gender.B.nz
NZ.conj$gender.B.nz.factor <- gender.B.nz.factor
NZ.conj$know.lgbt <- know.lgbt
NZ.conj$lgbt.self <- lgbt.self
NZ.conj$religious.yes <- religious.yes
NZ.conj$religiosity <- religiosity
NZ.conj$voted.yes <- voted.yes

NZ.conj$cm.left.leaning <- cm.left.leaning
NZ.conj$cm.pref.neighb <- cm.pref.neighb
NZ.conj$cm.electorab <- cm.electorab


#Subset of dataset
NZ.conj.labour <- subset(NZ.conj, party.choice.nz == "Labour")
NZ.conj.national <- subset(NZ.conj, party.choice.nz == "National")
NZ.conj.conserv <- subset(NZ.conj, soc.conserv.nz == 8 | soc.conserv.nz == 9 | soc.conserv.nz == 10)
NZ.conj.liber <- subset(NZ.conj, soc.conserv.nz == 0 |soc.conserv.nz == 1 | soc.conserv.nz == 2)
NZ.conj.women <- subset(NZ.conj, female.nz == 1)
NZ.conj.men <- subset(NZ.conj, female.nz == 0)
NZ.conj.knowlgbt <- subset(NZ.conj, know.lgbt == 1)
NZ.conj.notknowlgbt <- subset(NZ.conj, know.lgbt == 0)
NZ.conj.noreligserv <- subset(NZ.conj, religiosity == 0)
NZ.conj.religserv <- subset(NZ.conj, religiosity == 3 | religiosity == 4)
NZ.conj.young <- subset(NZ.conj, age.nz <= 35)
NZ.conj.old <- subset(NZ.conj, age.nz >= 60)
NZ.conj.LGBT.MP <- subset(NZ.conj, LGBT.all == 1)
NZ.conj.NO.LGBT.MP <- subset(NZ.conj, LGBT.all == 0)
NZ.conj.LGBT.MP.now <- subset(NZ.conj, LGBT.now == 1)
NZ.conj.NO.LGBT.MP.now <- subset(NZ.conj, LGBT.now == 0)
NZ.conj.lgbtvoters <- subset(NZ.conj, lgbt.self == 1)
NZ.conj.straightvoters <- subset(NZ.conj, lgbt.self == 0)



### Analysis - CONJOINT LGBT CANDIDATES


### Using [cregg] ###

#Rename conjoint attribute variables to have clearer plots and attach them to dataset 

Sexual_Orientation <- sexorient.B.nz.factor
Gender_Identity <- gender.B.nz.factor
Education <- educ.B.nz.factor
Race <- race.B.nz.factor
Political_Experience <- polexp.B.nz.factor
Religion <- religion.B.nz.factor
Age <- age.B.nz.factor
Health <- health.B.nz.factor

NZ.conj$Sexual_Orientation <- Sexual_Orientation
NZ.conj$Gender_Identity <- Gender_Identity
NZ.conj$Education <- Education
NZ.conj$Race <- Race
NZ.conj$Political_Experience <- Political_Experience
NZ.conj$Religion <- Religion
NZ.conj$Age <- Age
NZ.conj$Health <- Health

Choice <- dummy.choice.B.nz
NZ.conj$Choice <- Choice


#Create dummy variables for subgroup analysis and attach them to dataset

Labour = rep(NA,n) #vs. National
Labour[party.choice.nz=="Labour"]=1
Labour[party.choice.nz=="National"]=0 
table(Labour)
NZ.conj$Labour <- Labour

Progressive = rep(NA,n)
Progressive[soc.conserv.nz==0 | soc.conserv.nz==1 | soc.conserv.nz==2]=1
Progressive[soc.conserv.nz==8 | soc.conserv.nz==9 | soc.conserv.nz==10]=0
table(Progressive)
NZ.conj$Progressive <- Progressive

LGBT.friends <-know.lgbt
table(LGBT.friends)
NZ.conj$LGBT.friends <- LGBT.friends

Women<-female.nz
table(Women)
NZ.conj$Women <- Women

Young = rep(NA,n)
Young[age.nz<="35"]=1
Young[age.nz>="60"]=0
table(Young)
NZ.conj$Young <- Young

Non.religious = rep(NA,n)
Non.religious[religiosity==0]=1
Non.religious[religiosity==3 | religiosity==4]=0
table(Non.religious)
NZ.conj$Non.religious <- Non.religious


## Subset of dataset 

NZ.conj.labour <- subset(NZ.conj, party.choice.nz == "Labour")
NZ.conj.national <- subset(NZ.conj, party.choice.nz=="National")

table(party.choice.nz)



### ACME ###

#ACME - General
amce(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                 Health + Age + Religion + Political_Experience + 
                 Race + Education, id = ~ ResponseId)

amce.1 <- amce(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                 Health + Age + Religion + Political_Experience + 
                 Race + Education, id = ~ ResponseId)

plot(amce.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())


### MM ###

mm.1 <- mm(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
             Health + Age + Religion + Political_Experience + 
             Race + Education, id = ~ ResponseId, h0 = 0.5)

plot(mm.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) + geom_vline(xintercept = 0.50, colour="gray33")



## SUBGROUP ANALYSIS


# AMCE by groups (just change dataset in command)

#LGBT vs. straight voters

amce.lgbtvoters <- amce(NZ.conj.lgbtvoters, Choice ~ Sexual_Orientation + Gender_Identity +  
                 Health + Age + Religion + Political_Experience + 
                 Race + Education, id = ~ ResponseId)

amce.straightvoters <- amce(NZ.conj.straightvoters, Choice ~ Sexual_Orientation + Gender_Identity +  
                          Health + Age + Religion + Political_Experience + 
                          Race + Education, id = ~ ResponseId)



# MMs split by groups

#MM by gender
stacked.gender <- cj(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                Health + Age + Religion + Political_Experience + 
                Race + Education, id = ~ ResponseId,
              estimate = "mm", by = ~ Women)

#MM by party
stacked.gender <- cj(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                       Health + Age + Religion + Political_Experience + 
                       Race + Education, id = ~ ResponseId,
                     estimate = "mm", by = ~ Labour)

#MM by political ideology
stacked.gender <- cj(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                       Health + Age + Religion + Political_Experience + 
                       Race + Education, id = ~ ResponseId,
                     estimate = "mm", by = ~ Progressive)

#MM by having LGBT friends
stacked.gender <- cj(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                       Health + Age + Religion + Political_Experience + 
                       Race + Education, id = ~ ResponseId,
                     estimate = "mm", by = ~ LGBT.friends)

#MM by age
stacked.gender <- cj(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                       Health + Age + Religion + Political_Experience + 
                       Race + Education, id = ~ ResponseId,
                     estimate = "mm", by = ~ Young)

#MM by religiosity
stacked.gender <- cj(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                       Health + Age + Religion + Political_Experience + 
                       Race + Education, id = ~ ResponseId,
                     estimate = "mm", by = ~ Non.religious)

##MM by sexual orientation
stacked.sexor <- cj(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                       Health + Age + Religion + Political_Experience + 
                       Race + Education, id = ~ ResponseId,
                     estimate = "mm", by = ~ lgbt.self)


### SUBGROUP DIFFERENCES ###

### ACME differences
acme.diff.1 <- amce_diffs(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, ~Women , id = ~ ResponseId)

plot(acme.diff.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())


### MM differences

#mm differences by gender
mm.diff.gender <- mm_diffs(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                             Health + Age + Religion + Political_Experience + 
                             Race + Education, 
                           ~ Women, id = ~ ResponseId)
plot(mm.diff.gender, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by party
mm.diff.party <- mm_diffs(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                            Health + Age + Religion + Political_Experience + 
                            Race + Education, 
                          ~ Labour, id = ~ ResponseId)
plot(mm.diff.party, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by political ideology
mm.diff.polideol <- mm_diffs(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                               Health + Age + Religion + Political_Experience + 
                               Race + Education, 
                             ~ Progressive, id = ~ ResponseId)
plot(mm.diff.polideol, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by having LGBT friends
mm.diff.lgbtfriends <- mm_diffs(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                                  Health + Age + Religion + Political_Experience + 
                                  Race + Education, 
                                ~ LGBT.friends, id = ~ ResponseId)
plot(mm.diff.lgbtfriends, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by age
mm.diff.age <- mm_diffs(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                          Health + Age + Religion + Political_Experience + 
                          Race + Education, 
                        ~ Young, id = ~ ResponseId)
plot(mm.diff.age, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by religiosity
mm.diff.religiosity <- mm_diffs(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                                  Health + Age + Religion + Political_Experience + 
                                  Race + Education, 
                                ~ Non.religious, id = ~ ResponseId)
plot(mm.diff.religiosity, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank())

#mm differences by sexual orientation
mm.diff.sexor <- mm_diffs(NZ.conj, Choice ~ Sexual_Orientation + Gender_Identity +  
                                  Health + Age + Religion + Political_Experience + 
                                  Race + Education, 
                                ~ lgbt.self, id = ~ ResponseId)





## Analysis for causal mechanism



#Causal mediation

library(mediation)

med.fit <- lm(cm.pref.neighb ~ gender.B.nz.factor + sexorient.B.nz.factor + health.B.nz.factor +
                 age.B.nz.factor + religion.B.nz.factor + polexp.B.nz.factor + 
                 race.B.nz.factor + educ.B.nz.factor +
                 age.nz + educ.nz + female.nz + know.lgbt + party.choice.nz +
                 religiosity + income.h.nz, data = NZ.conj)
out.fit <- lm(dummy.choice.B.nz ~ cm.pref.neighb + gender.B.nz.factor + sexorient.B.nz.factor + health.B.nz.factor +
                age.B.nz.factor + religion.B.nz.factor + polexp.B.nz.factor + 
                race.B.nz.factor + educ.B.nz.factor +
                age.nz + educ.nz + female.nz + know.lgbt + party.choice.nz +
                religiosity + income.h.nz, data = NZ.conj) 

med.out <- mediate(med.fit, out.fit, treat = "sexorient.B.nz.factor", 
                   mediator = "cm.pref.neighb", 
                   sims = 1000,
                   dropobs = TRUE,
                   robustSE = TRUE)
summary(med.out)
plot(med.out)


#Sensitivity analysis
sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out)


## CAUSAL MECHANISMS with mechanisms as OUTCOME variables

#Electability
amce.cm.1 <- amce(NZ.conj, cm.electorab ~ Sexual_Orientation + Gender_Identity +  
                 Health + Age + Religion + Political_Experience + 
                 Race + Education, id = ~ ResponseId)

plot(amce.cm.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) +
  xlim(-0.15, 0.16)


#Prejudice
amce.cm.2 <- amce(NZ.conj, cm.pref.neighb ~ Sexual_Orientation + Gender_Identity +  
                 Health + Age + Religion + Political_Experience + 
                 Race + Education, id = ~ ResponseId)

plot(amce.cm.2, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) +
  xlim(-0.15, 0.16)


#Ideology - left
amce.cm.3 <- amce(NZ.conj, cm.left.leaning ~ Sexual_Orientation + Gender_Identity +  
                 Health + Age + Religion + Political_Experience + 
                 Race + Education, id = ~ ResponseId)

plot(amce.cm.3) + theme(legend.title = element_blank())

plot(amce.cm.3, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme(legend.title = element_blank()) +
  theme(strip.text = element_blank()) +
  xlim(-0.15, 0.16)




## INTERACTION EFFECTS: FINDIT PACKAGE

#Subset so you don't have missing data
myvars2 <- c("dummy.choice.B.nz", "religion.B.nz.factor", "gender.B.nz.factor", "sexorient.B.nz.factor",
             "health.B.nz.factor", "educ.B.nz.factor", "age.B.nz.factor", "religion.B.nz.factor",
             "polexp.B.nz.factor", "race.B.nz.factor", 
             "ResponseId")
NZ.conj.LGBT.inter <- NZ.conj[myvars2]
NZ.conj.LGBT.inter <- na.omit(NZ.conj.LGBT.inter)

## Specify the order of each factor
NZ.conj.LGBT.inter$gender.B.nz.factor<- factor(NZ.conj.LGBT.inter$gender.B.nz.factor,ordered=TRUE,
                                               levels=c("Woman", "Transgender", "Man"))
NZ.conj.LGBT.inter$sexorient.B.nz.factor<- factor(NZ.conj.LGBT.inter$sexorient.B.nz.factor,ordered=TRUE,
                                                  levels=c("Straight", "Gay"))
NZ.conj.LGBT.inter$race.B.nz.factor<- factor(NZ.conj.LGBT.inter$race.B.nz.factor,ordered=TRUE,
                                             levels=c("Maori-Pacific", "White"))

#Calculate AMIE 

#Gender and sexual orientation 
fit.1 <- CausalANOVA(formula=dummy.choice.B.nz ~ gender.B.nz.factor + sexorient.B.nz.factor + 
                       health.B.nz.factor +
                       age.B.nz.factor + religion.B.nz.factor + polexp.B.nz.factor + 
                       race.B.nz.factor +
                       educ.B.nz.factor,
                     int2.formula = ~ gender.B.nz.factor:sexorient.B.nz.factor,
                     data=NZ.conj.LGBT.inter, 
                     cluster=NZ.conj.LGBT.inter$ResponseId, nway=2)
summary(fit.1)

#Sexual orientation and race
fit.2 <- CausalANOVA(formula=dummy.choice.B.nz ~ gender.B.nz.factor + sexorient.B.nz.factor + 
                       health.B.nz.factor +
                       age.B.nz.factor + religion.B.nz.factor + polexp.B.nz.factor + 
                       race.B.nz.factor +
                       educ.B.nz.factor,
                     int2.formula = ~ sexorient.B.nz.factor:race.B.nz.factor,
                     data=NZ.conj.LGBT.inter, 
                     cluster=NZ.conj.LGBT.inter$ResponseId, nway=2)
summary(fit.2)


#Gender and race
NZ.conj.LGBT.inter$gender.B.nz.factor<- factor(NZ.conj.LGBT.inter$gender.B.nz.factor,ordered=TRUE,
                                               levels=c("Woman", "Man", "Transgender"))
fit.3 <- CausalANOVA(formula=dummy.choice.B.nz ~ gender.B.nz.factor + sexorient.B.nz.factor + 
                       health.B.nz.factor +
                       age.B.nz.factor + religion.B.nz.factor + polexp.B.nz.factor + 
                       race.B.nz.factor +
                       educ.B.nz.factor,
                     int2.formula = ~ gender.B.nz.factor:race.B.nz.factor,
                     data=NZ.conj.LGBT.inter, 
                     cluster=NZ.conj.LGBT.inter$ResponseId, nway=2)
summary(fit.3)


#Health and race
NZ.conj.LGBT.inter$health.B.nz.factor<- factor(NZ.conj.LGBT.inter$health.B.nz.factor,ordered=TRUE,
                                               levels=c("HIV+ since birth", "Weelchair since birth", 
                                                        "Healthy", "Overweight, diabetes", "HIV+"))
fit.4 <- CausalANOVA(formula=dummy.choice.B.nz ~ gender.B.nz.factor + sexorient.B.nz.factor + 
                       health.B.nz.factor +
                       age.B.nz.factor + religion.B.nz.factor + polexp.B.nz.factor + 
                       race.B.nz.factor +
                       educ.B.nz.factor,
                     int2.formula = ~ health.B.nz.factor:race.B.nz.factor,
                     data=NZ.conj.LGBT.inter, 
                     cluster=NZ.conj.LGBT.inter$ResponseId, nway=2)
summary(fit.4)


